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The ground state parameters of the two-dimensional S = 1/2 antiferromagnetic Heisenberg 
model are calculated using the Stochastic Series Expansion quantum Monte Carlo method for L x L 
lattices with L up to 16. The finite-size results for the energy E, the sublattice magnetization 
M, the long- wavelength susceptibility X-i-il = 27r/L), and the spin stiffness ps, are extrapolated to 
the thermodynamic limit using fits to polynomials in 1/L, constrained by scaling forms previously 
obtained from renormalization group calculations for the nonlinear a model and chiral perturbation 
theory. The results are fully consistent with the predicted leading finite-size corrections and are of 
sufficient accuracy for extracting also subleading terms. The subleading energy correction (~ i/L'^) 
agrees with chiral perturbation theory to within a statistical error of a few percent, thus providing 
the first numerical confirmation of the finite-size scaling forms to this order. The extrapolated 
ground state energy per spin, E = —0.669437(5), is the most accurate estimate reported to date. 
The most accurate Green's function Monte Carlo (GFMC) result is slightly higher than this value, 
most likely due to a small systematic error originating from "population control" bias in GFMC. The 
other extrapolated parameters are AI = 0.3070(3), ps = 0.175(2), x± = 0.0625(9), and the spinwave 
velocity c = 1.673(7). The statistical errors are comparable with those of the best previous estimates, 
obtained by fitting loop algorithm quantum Monte Carlo data to finite-temperature scaling forms. 
Both AI and ps obtained from the finite-T data are, however, a few error bars higher than the 
present estimates. It is argued that the T = extrapolations performed here are less sensitive to 
effects of neglected higher-order corrections and therefore should be more reliable. 



I. INTRODUCTION 

In the nonlinear a model description of the two-dimensional (2D) Heisenberg modelj3 the low-energy and low- 
temperature properties of the system are completely determined by three ground state parameters; the sublattice 
magnetization M , the spin stiffness constant p^, and the spinwave velocity c. Their values are not given by the theory, 
however, but have to be determined starting from the microscopic Hamiltonian. A large number of calculations of 
the ground state parameters have been_carried out. The antiferromagnetically ordered ground state, which has been 
established rigorously only for S > W2,B was first convincingly confirmed also for 5* = 1/2 in a quantum Monte Carlo 
(QMC) study by Reger and Young.S The sublattice magnetizatian obtained this way, M 0.30 (in units where the 
Neel state has M = 1/2), also indicated that spinwave theoryafl Sfves a surprisingly good quantitative description 
of the ground state. The same conclusion was reached by Singh,El who carried out a series expansion around the 
Ising limit, and found M « 0.30, ps ~ 0.18J, and c w 1.7J (^|4s the nearest-neighbor exchange coupling), all in good 
agreement with spinwave theory including the 1/S correctjLj)^g|.l3 Subsequent higher-order spin W|aj|ja-| ^( juj|^jj^^i.|h|^ 

ojje^confirmed and improved on the 
indicate that the true values of 
the ground state parameters deviate from their l/S*^ spinwave values by only 1-2% or less. 

For most practical purposes (such as extracting J for a system from experimental data) , the ground state parameters 
of the 2D Heisenberg model are now known to quite sufficient accuracy. However, there are still reasons to go to 
even higher precision. One is that the model is one of the basic "prototypic" many-body models in condensed matter 
physics. It has become a testing ground for various analytical and numerical methods for strongly correlated systems, 
thus making it important to accurately establish its properties. Another reason is the very detailed predictions 
that haySj. r i j i su lted from field theoretical studies, such as renormalization group calculations for the nonlinear a 



model,Elcil'E3c3 and chiral perturbation theory.cZl Apart from giving the low-energy projpe.rties in the thermodynamic 
limit, these theories also predict the system size dependence of various quantities This is important from 
the standpoint of numerical calculations such as exact diagonalization and QMC, which are necessarily restricted 
to relatively small lattices. Finite-size scaling approaches have been very successful in the study of ID quantum 
spin systems, having convincingly confirmed various predictions from bosonization and conformal field theory. For 
example, critical exponents aad logarithmic corrections have been extracted from the size dependence of ground state 
energies and finite-size gapsJHa and from correlation functions.^ With the concrete predictions now available, similar 



that the 1/5^ correctioBS.tQ. M, ps and c indeed are smallE'Ercl Several other QMC simulations 
exact diagonalizationSjcScilO as well as series expansions to higher ordera£3_hai£e.confirme 
accuracy of the above estimates. The presently most accurate calculationgl3ll3'ESt3 
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studies show great promise for testing theories also in 2D. For the standard Heisenberg model, finite-size scaling has 
been used esteft^yelw-^ffld successfully in extrapolating, e.g., the sublattice magnetization for small lattices to infinite 
system sizeP't2Hlilll3't3EJ but only a few studies haveusctikci jb^en accurate enough for reliably addressing the validity 



of the theoretical predictions for the size corrections. L-ii—ii—ii—i ■— ■ 

In one dimension, exact diagonalization, and more recently the density matrix renormalization group mcthod,E2l 
enable highly ajccucate calculations for systems sufficiently large to approach the limit where the asymptotic scaling 
forms are validH3'E3 Calculations with these methods in two dimensions cannot reach linear dimensions large enough 
to verify the details of the predicted scaling forms, however. Some of the expected ieading finite-size behavior has 
been seen in exact diagonalization studies including systems with up to 6 x 6 spinspJ^Ell but constants extracted from 
the size dependence are typically not consistent with other calculations. For example, c extracted from the scaling of 
E deviates by 15% from other estimates. E3 There are hence clear indications that these small systems are not yet in 
the regime where only the dominant corrections are important. 

QMC can reach significantly larger lattices at the cost of statistical errors which are often relatively large, making 
it difficult to accurately extract the scaling behavior. Runge carried out Green's function Monte Carlo (GFMC) 
simulations of L x L Heisenberg systems with L up to 16, and fo und a reasonable consistency with the leading 
T — size dependence of the energy and the sublattice magnetization.llaO He also noted the presence of a subleading 
correction to the energyjlj but the accuracy of the GFMC data was not high enough to extract its order, and 
furthermore c extracted from the leading correction was sensitive to the subleading one. The_extrapolated ground 
state energy obtained in this study is nevertheless the most accurate estimate reported so far.EEl 

Chiral perturbation theory has recently enabled calculations of scaling forms for finite size and finite temperature fpj: 
various quantities Ej Such forms have been used in combination with QMC data in recent work by Wiese and Ying,E£l 
and Beard and Wiese.E3 Their calculations employed, respectively, methods based on the "loop-cluster algO|rithm" 
suggested by Evertz et alS^, and a continuous-time variant of that method developed by Beard and Wiese.t2l These 
algorithms are based on global fL ia s of loops of spins, and overcome the problems with long autocorrelation times 
typical of standard Suzuki- TrotterSO or worldlineE^ QMC methods (the continuous-time approach furthermore avoids 
the systematic discretization error of the Trotter approximation). Considerably more accurate finite-T data could 
therefore be generated, and the leading-order scaling forms of chiral perturbation theory were convincingly verified, 
both in the "cubic" regime T/c « 1/L (Ref. |l^) and the "cylindrical" regime T/c <C 1/i (Ref. |l^). The extrapolated 
Af , ps and c are the most accurate reported to date, although there are some minor discrepancies between the two 
results for M (on the border line of what could be expected within statistical errors alone). 

In this paper, a finitfigsize scaling study of T = data is reported. Using the the Stochastic Series Expansion 
(SSE) QMC algorithmJ^'EEl energy results of unprecedented accuracy are obtained for L x L lattices with L up to 
16. The relative statistical errors are as low as « 10~^. limploying a recently suggested data analysis scheme which 
takes into account covariance among calculated quantiticSjEZI very accurate results for the sublattice magnetization are 
also obtained. Furthermore, the spin stiffness and the long-wavelength susceptibility x(9 — 27t/L) are also calculated 
directly in the simulations. 

Assuming for the size depjsndences polynomials in 1/L, constrained by scaling forms for E and AI predicted from 
chiral perturbation theory,E3 all the computed quantities are included in a coupled fit- The quality of the QMC 
data for E and AI is high enough that size corrections beyond the subleading terms have to be included. The leading- 
order corrections are fully consistent with the predictions. From a careful statistical analysis of the fits, bounds for 
the subleading terms are estimated. The subleading energy correction is found to agree with the prediction of chiral 
perturbation theory to within a statistical error of 5% (the subleading correction for M is also estimated, but has not 
yet been calculated analytically). This is the first numerical confirmation of chiral perturbation theory to subleading 
order. 

The extrapolated ground state energy, E = —0.669437(5), is the-jnost accurate estimate reported to date, with a 
statistical error six times smaller than the GFMC result by Runge,E£l and is slightly lower than his result. Comparing 
the finite-size data of the two calculations, a clear tendency to over-estimation of the energy is seen in the GFMC 
results. This is likely due to |£ubias originating from "population control" in GFMC (a small effect of this nature was 
in fact anticipated by RungetS). 

The results for the sublattice magnetization, M — 0.3070(3), and the spin-stiffness, ps — 0.175(2), are both slightly 
lower than the estimates from the finite-T scaling by Beard and WiesetJ [M — 0.3083(2) and ps — 0.185(2)]. Although 
it is at this point difficult to definitely conclude which calculation is more reliable, it can again be noted that the high 
accuracy of the QMC data for E and M used in the fits carped out in this paper necessitates the inclusion of size 
corrections beyond the orders considered by Beard and Wiese.Ej Hence, any remaining effects of neglected corrections 
of even higher order should be smaller. Quantitative estimates of such effects on the extrapolations performed here 
support that they are smaller than the statistical errors indicated above. For the spinwave velocity the two calculations 
agree, the result obtained here being c — 1.673(7) and the value reported in Ref. 19 being c = 1.68(1). 

The outline of the rest of the paper is the following. In Sec. II the SSE algorithm and the covariance error reduction 
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scheme are outlined. The absence of systematic errors are demonstrated in comparisons with exact results for 4 x 4 
and 6x6 lattices. The fitting procedures and the results of these are discussed in Sec. III. The study is summarized 
in Sec. IV. Some other problems where the methods applied here should be useful are also mentioned. 



II. NUMERICAL METHODS 

The standard 2D Heisenberg model is defined by the Hamiltonian 

H = JJ2^^■^J^ (J>0), (1) 

where Si is a spin-1/2 operator at site i on a square lattice with N — L x L sites, and denotes a pair of 

nearest-neighbor sites. Below, in II-A, the SSE approach to QMC simulation of this model is outlined. More details 
of the ^gi^^ are discussed in Refs. 36 and |3^. The SSE rr^hod has recently been applied to a variety of spin 



models,t as well as ID Hubbard- type electronic models O 

It was recently noted that correlations between measurements of diffeijeait observables can be used to significantly 
increase the accuracy of certain quantities calculated in SSE simulations .l3 This covariance scheme for analyzing the 
data is of crucial importance in the present work, and therefore this method is also described below in Sec. II-B. The 
high accuracy of the procedures is demonstrated by comparing results for 4x4 and 6x6 lattices with the exact 
diagonalization data available for these systems. 

A. Stochastic Series Expansion 

Based on the exact powcE-secies expansion of e-^", the SSE metho dii can be considered a generalization of 
Handscomb's QMC scheme. ejO It is the first "exact" method proposed for QMC simulations of general lattice 
Hamiltonians at finite-temperature (with the usual caveat of being in practice restricted to models for which the sign 
problem can be avoided). It-4s hence not based on a controlled approximation, such as the Trotter formula used in 
standard worldline methods,Ell and directly gives results accurate to within statistical errors. Despite being formulated 
at finite T, temperatures low enough for studying the ground state can easily be reached for moderate-size lattices. 

As in Handscomb's method for the 5 = 1/2 antiferromagnet,c3 the SSE approach for this model starts from the 
Hamiltonian written as {J — 1) 



H 



-^[Hl,fa-i/2,b] + y, (2) 
fc=l 



where 6 is a bond connecting a pair of nearest-neighbor sites (i(b),j{b)), and the operators Hi^b and J?2,ti are defined 
as 

= 2[| — S'j^(j)S'|(b)] (3a) 

H2.b = S^(^b)^J'{b) + ^i{b)^t{b)- 

An exact expression for an operator expectation value 

(A) = ^TT{Ae'f'"}, Z^Tr{e-^"}, (4) 

at inverse temperature /3 = J/T, is obtained by Taylor expanding e~^^ and writi^ the traces as sums over diagonal 
matrix elements in the basis {|a) = \Si, . . . , S^)}. The partition function is thcnl 



a 



i=l 



(5) 



a n=0 S„ 

where Sn is a sequence of index pairs defining the operator string Y['i=i ^ai,bi, 

5„ = [ai,6i][a2,62]...[a„,6„], a, G {1, 2}, 6, e {1, . . . , 2A^}, (6) 
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and 722 denotes the totajLuumber of index pairs (operators) [oi, bi] with = 2 (n = rii + 712). Eq. (||) deviates from 
Handscomb's niethod,OE3 which rehes on exact evaluation of the traces of the operator sequences and therefore is 
hmited to models for which this is possible. The Heisenberg model considered here is such a model (one of the very 
few), but the more general SSE approach of further expanding over a set of basis states is preferable also in this case, 
for reasons that will be discussed below. 

The objective now is to develop a scheme for importance sampling of the terms in the partition function (^). A 
term, or configuration, (a, S'„) is specified by a basis state \a) and an operator sequence Sn- The operators Hi^t, and 
H2,b can act only on states where the spins at sites i{b) and j{b) are antiparallel. The diagonal Hi t leaves such a 
state unchanged, whereas the off-diagonal i?2,b flips the spin pair. Defining a propagated state 

p 

Hp)) = Y[^-^M, KO)) = |a), (7) 

a configuration {a,Sn) must clearly satisfy the periodicity condition |a(n)) = |q^(0)) in order to contribute to the 
partition function. For a lattice with L x L sites and L even, this implies that the total number 712 of the off-diagonal 
operators must be even, and hence that all terms in (^) are positive and can be used as relative probabilities in a 
Monte Carlo importance sampling procedure (this is true for any non-frustrated system). 

For a finite system at finite /3, the powers n contributing significantly to the partition function are restricted to within 
a well defined regime, and the sampling space is therefore finite in practice. In order to construct an efficient updating 
scheme for the index sequence it is useful to explicitly truncate the Taylor expansion at some self-copaistently chosen 
upper bound n ~ I, high enough to cause only an exponentially small, completely negligible error.ta One can then 
define a sampling space where the length of the sequence is fixed, by inserting a number l — n of unit operators, denoted 
Hofi, in the operator strings. The terms in the partition function (|^) are divided by (^V, in order to compensate for 
the number of different ways of inserting the unit operators. The summation over n in (^ is then implicitly included 
in the summation over all sequences Si of length I, with [a^, bi] = [0, 0] as an allowed operator. Denoting by W{a, Si) 
the weight of a configuration {a, Si), the partition fimction is then 

a Si 

Since all non-zero matrix elements in (^ equal one, the weight is (when non- vanishing) 

,n«,5,)^(f)"<^I^^. (9) 

where n still is the expansion power of the term, i.e., the number of non-[0, 0] operators in Si. 

The following is only a brief outline of the actual sampling scheme. More details can be found in Rcfs. ^ and |3|. 
During the simulation. Si and one of the states \a{p)) are stored. Other propagated states are generated as needed. 
The simulation is started with a randomly generated state |q;(0)), with an index sequence Si containing only [0,0] 
operators (unit operators), and with some arbitrary (small) I. The truncation I is adjusted as the simulation proceeds, 
as will be discussed further below. 

With the fixed-length scheme, all updates of the operator sequence can be formulated in terms of substitutions of 
one or several operators. The simplest involves a diagonal operator at a single position; [0,0] ^ [li^]- This update 
can be carried out consecutively at all positions p for which Up G {0, 1}. In the — > direction, the bond index b is 
chosen at random and the update is rejiectcd if the spins connected by b are parallel in the current state \a{p — 1)). 
The Metropolis acceptance probabilitietO required to satisfy detailed balance are obtained from (^, where the power 
n in is changed by ±1. Updates involving the off-diagonal operators [2, b] are carried out with n fixed. The simplest is 
of the type [1, b][l, b] [2, 6] [2, b], involving two operators acting on the same bond. These two sequence updates can 
generate all configurations with spin flips on retracing paths on the lattice, and are the only ones required for a ID 
system with open boundary conditions. For a 2D system, configurations associated with spin flips around any closed 
loop are possible, and an additional type of update is required. It is sufficient to consider substitutions on a plaquette, 
of the type [2, bi] [2, 62] ^ [2, 63] [2, 64], where 61, . . . , 64 is a permutation of the four bonds of a plaquette. For systems 
with periodic boundary conditions, updates involving cyclic spin flips on loops wrapping around the whole system 
are required (sampling of different winding number sectors), and cannot be accomplished by the above local sequence 
alterations. For the square lattice considered here, the winding number can be changed by substituting L/2 operators 
according to [2, 61] . . . [2, 6^/2] <-» [2, 6^/2+1] . . . [2, 6l], where the set of bonds &i, 6l is a permutation of bonds 
forming a closed ring around the system in the x- or j/-direction. 
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Updating the operator sequence with the four types of operator substitutions described above suffices for generating 
all possible configurations within a sector of fixed magnetization, = ^iLi ■ the grand canonical ensemble, 
global spin flips changing the magnetization are also required— Here T — )■ will be considered (i.e., T is much lower 
than the finite-size gap), and since the ground state is a singletcZI the canonical ensemble with = is appropriate. It 
can be noted that in Handscomb's method the sampling is (in principle) automatically over all magnetization sectors, 
and therefore a restriction to, e.g., = is not possible. In practice, this causes problems at low temperatures, and 
Handscomb's method has thereface been used for the antiferromagnetic Heisenberg model mostly at relatively high 
temperatures (T/J > 0.4 in 2D).LJ The SSE method with the restriction = can be used at arbitrarily low T. 

In order to determine a sufficiently high truncation of the expansion, the fluctuating power n is monitored during 
the equilibration part of the simulation. If n exceeds some threshold I — Ai/2, the cut-off is increased, I ^ I + Ai, by 
inserting additional -ffo,o operators at random positions. In practice A; « I/IO leads to a rapid saturation of ^ at a 
value sufficient to cause no detectable truncation errors. The growth of I during equilibration is illustrated for a 4 x 4 
system in Fig. ^ The distribution of n during a subsequent simulation is shown in Fig. |2[ and clearly demonstrates 
that the truncation of the expansion is no approximation in practice. 

A Monte Carlo step (MC step) is defined as a series of the single-(diagonal)operator substitutions attempted 
consecutively at each position in Si (where possible), followed by a series of off-diagonal updates carried out on each 
bond, plaquette, and ring. Due to the locality of the, constraints in these updates, the number of operations (the CPU 
time) per MC step scales linearly with N and However, the acceptance rate for the "ring update" that changes 

the winding number decreases rapidly with increasing system size. It is therefore sometimes useful to increase the 
number of attempted ring updates with the system size, which then leads to a faster growth of CPU time with N. 
The acceptance rate of the ring update currently used becomes too low for L > 16, and simulations of larger aistems 
therefore in practice have to be restricted to the sector with zero winding number. It has recently been notedCj that 
in fact exact results are obtained as T ^ even for simulations restricted this way. However, compared to simulations 
with fluctuating winding numbers, lower temperatures are required for the system observables to saturate at their 
ground state values. e2I Here only systems with L < 16 are considered, and the update changing the winding number 
is always included. 

Measurements of physical observables are carried out using the index sequences Sn obtained by omitting the [0, 0] 
operators in the generated Si. These are then, of course, distributed according to the weight function corresponding 
to Eq. d). 

One can show that .the internal energy per spin is simply given by the average of n [with the constant term in 
Eq. (I) neglected] :m 

This expression also shows that the average power, and hence the sequence length I, scales as f3N at low temperatures. 
A spin-spin correlation function, 

C(z,j)=C(r,-r,)==(5f5|), (11) 
is obtained averaging the correlations in the propagated states \a{p)) defined in Eq. (^. Further defining 

Snp] = {a{p)\Sna{p)), (12) 

the correlation function is given byll 

c^(^j)-(;^E^fW5|b]). (13) 

The corresponding static susceptibility, 

X{t,j)= I dT{S:{T)S^{0)), (14) 



p=0 







involves correlations between all the propagated states :cj 
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OfT-diagonal correlation functions can be easily calculated for operators that can be expressed in terms of the spin- 
flipping operators H^fi^ each of which is a sum of two terms; H'^ = Sff^ij^Sjfj^-^ and = S~^^^-^S'^(j^y The spin stiffness 
constant involves a static susceptibility defined in terms of these operators. 

Although the simulation scheme is formulated with H2.b — + , one can still access the terms individually 
since only one of them can propagate a given state. One can show that an equal-time correlation function, 

F,,,ib,b') = {HSHS,'), (16) 

is given by@ 

F„.,ib,b')^(^^j^Niba;b'a')), (17) 



where N(ba; b'a') is the number of times the operators and H^, appear next to each other in Sn, in the given 
order. The corresponding static susceptibility, 

13 

Xaa'{b,b')= /dr(K(r)K''(0)), (18) 



is given by the remarkably simple formulate 

Xaa' [b, b') = 4 {N{ba)N{b'a') - 5w5..'N{ba)) (19) 

where N{ba) is the total number of operators in Sn- 

Now a direct estimator for the spin stiffness can be constructed. The stiffness, ps, is defined as the second derivative 
of the ground state energy with respect to a twist $ in the boundary condition, around an axis perpendicular to the 
direction of the broken symmetry. For a finite lattice, where the symmetry is not broken, a factor 3/2 has to be 
included in order to account for rotational averaging. Distributing the twist equally over all interacting spin pairs 
{i,j)x in the a;-direction, the finite-size definition for ps is hence 



3 1 d^EnicI)) 

Pi - 



2L2 



, (20) 

0=0 



where = ^/L. An expression which is only dependent on the ground state at = is obtained by expanding the 
Hamiltonian to second order in (p. The Hamiltonian in the presence of the twist is 



77(0)= S..i?(0)S,+ S.-S,, (21) 



where R{(j)) is the rotation matrix 



cos {4>) sin [(j)) 

R{(l)) =1 - sin {(j)) cos (0) I . (22) 
1 

Expanding to second order in results in 

i?(0) - H(0) = E WiSfS^^ + Sfsy) + id,{StSj - Srs+)] . (23) 

The first term is proportional to H{Q) (for the rotationally invariant case considered here). The expectation value of 
the second term vanishes, but it gives a contribution quadratic in in second order perturbation theory. Defining the 
spin current operator 

j^ = l E (^^^i - ^r^;), (24) 

and the current-current correlation function at Matsubara frequency w„i — 27rmT, 
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= dTe~''^"^^js{T)jM), (25) 



the stiffness is given by 

Ps = -li^E + AsiO)], (26) 

wfiere E is tlie ground state energy per spin. 

Tfie QMC estimate for tlie energy is given by Eq. ([l0|). The current-current correlator = As(0) is a sum of 
integrals of the form (p^). Denoting by and the number in Sn of operators S^Sj' and S^Sj' with a 
bond in the x-direction, Eqs. ( |25| ) and ( p^ give 

Ps = |§((A^.+ -iV-)2), (27) 

i.e. the terms linear in iV+ and cancel. Defining the winding numbers Wx and Wy in the x and y direction: 

= (iV+ - N-)/L, {a = y), (28) 

the stiffness can also be written as 

p, = l(wl+wl)/p. (29) 

This definition is clearly valid only for a simulation that samples all winding number sectors. With a restriction to 
the subspace with Wx — Wy — 0, ps can be calculated using the long-wavelength limit of a current-current correlator 
involving a twist field with a spatial modulation.cS 

The above method of calculating the stiffness directly from the winding numbe|r fluctuations is clearly strongly 
related to methods used for the superfluid density in simulations of boson models.E3 



B. Error Reduction Using Covariance 



In Monte Carlo simulations, fluctuations (statistical errors) of different physical observables are often correlated with, 
each other. These covariance effects can sometimes be used to obtain improved estimators for certain quantities.Ell 
In some cases one may have exact knowledge of some quantity independently of the QMC calculation. If there are 
strong correlations between a known quantity A and some other, unknown quantity B, the accuracy of B can be 
improved via its covariance with the measured A, by calculating the average and statistical error under the condition 
that A equals its known value. In other cases, it may be possible to calculate a quantity in more than one way in 
the same simulation. If one of the estimates, Ai, is more accurate than the other, A2, a covariance between A2 and 
some other quantity B can again be used to improve the estimate of B. With the SSE method the internal energy 
of the rotationally invariant Heisenberg model can be calculated in two different ways: Ei from the average power 
of the series expansion according to Eq. (p^), and using the nearest-neighbor correlation function C(1,0) calculated 
according to Eq. (^3|); E2 — 6C(1,0). The manifestly rotationally invariant estimator Ei is signiflcantly less noisy 
than i?2. Results for quantities with fluctuations correlated to those of C(l, 0), such as C(r) with r > 1, can therefore 
be improved with the aid of Ei . 

For the purpose of accurately measurins-correlations between the fluctuations of two different quantities, the so 
called "bootstrap method" is a useful toolEl With the simulation data as usually divided into M "bin" averages, a 
"bootstrap sample" A^ is defined as an average over M randomly selected bins (i.e., the same number as the total 
number of bins, allowing, of course, multiple selections of the same bin). With r{i) denoting the i:th randomly chosen 
bin, 

1 



The statistical error can be calculated on the basis of A/^ bootstrap samples Aji^ , according ti 




, Mr 
i—l 
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where A is the regular average over all bins. Note that Eq. (|3|) lacks the factor (M/j — present in the conventional 
expression for the variance of the average calculated on the basis of Mn bins. The bootstrap method is in general 
more accurate (due to a better realization of a Gaussian distribution for the bootstrap samples), in particular if A 
is not measured directly in the simulation, but is some nonlinear function of measured quantities (in which case A 
should be calculated on the basis of bootstrap samples, not individual bins). Sets of bootstrap samples and 
{Br-} generated on the basis of the same randomly selected bins are well suited for evaluating correlations between 
the statistical fluctuations of A and B, and are used in the covariance error reduction scheme described next. 

Here this method will be illustrated using simulation results for the staggered structure factor S'(7r, tt) and the 
staggered susceptibility xi'^,'^)- These are defined according to 

'^(^''^) = ^E(-1)'^'""''"^"''C'(^,J) (32a) 

with C{i,j) and xihj) given by Eqs. ( p^ and ([T5|). The structure factor is of particular interest, since it defines the 
sublattice magnetization squared of a finite system. The fluctuations of S{tt, tt) are strongly correlated to those of 
C(1,0), and S'(7r, tt) can therefore be calculated to an accuracy significantly higher than if only the direct estimator 
( p^ ) is used. The susceptibility x(7r,7r) is only weakly correlated with C(1,0), however, and only a modest gain in 
accuracy can be achieved for this quantity. 

First some results for a 6 x 6 lattice are discussed. This is the largest system for which Lanczos results have been 
obtainedH Comparing with these exact results, the accuracy of the QMC technique and the covariance method can 
be rigorously checked. The temperature used in the simulation has to be low enough for the calculated quantities 
to have saturated at their ground state values. In order to check for temperature effects, several calculations were 
carried out. Results at inverse temperatures /3 — 2A and 48 are indistinguishable within error bars, indicating that 
these temperatures are sufficiently low for L — 6. The results presented below are for (3 — 48. The simulation was 
divided into bins of 5 x 10^ MC step each, and a total of 600 bins were generated. 

Fig. 1^ shows the covariance between the measured nearest-neighbor correlation function {E2) and ^(Tr, tt). The 
plot was generated on the basis of 2000 bootstrap samples. Strong linear correlations between the two quantities are 
evident. Hence further knowledge of E can improve the estimate of S. The conventional average and error of S'(7r, tt) 
is calculated on the basis of all the points i.e., the distribution obtained by projecting the points onto the 5'-axis. 
Having a better estimate Ei ± ai for E, an improved estimate of S can be calculated by weighting the points by a 
Gaussian centered at Ei and with a width equal to the error ai. In this case, the reduced statistical error is w 1/12 of 
the conventional error. Note that the conventional estimates of both S and E2 lie outside the exact results by « 1.5 
standard deviations (not an unlikely situation statistically). The improved estimate of S is nevertheless within one 
standard deviation of the exact result, reflecting this being the case for the more accurate energy estimate Ei used 
in the procedure. In fact, this correcting property of the covariance method can even eliminate certain systematic 
errors, such as those originating from finite-T effects in calculations aimed at ground state properties.E3 

Besides illustrating the use of the covariance method to reduce the statistical errors. Fig. |^ also clearly demonstrates 
to a high accuracy the absence of detectable systematic errors in the QMC data. This confirms that the SSE method 
indeed produces unbiased results. Table || summarizes the comparisons with the exact results for both 4x4 and 6x6 
lattices. 

As the system size increases, the fluctuations in S{'k, tt) as computed in the standard way increase, and accurate 
estimates become increasingly difficult to obtain. This is typical of algorithms utilizing local updates. The fluctuations 
in the energy per site as calculated from {n) actually decrease, however (due to self-averaging). Hence, the gain in 
accuracy achieved with the covariance effect increases with the system size. Fig. ^ shows L = 16 results for the 
staggered structure factor. For this system size the error in the energy estimate Ei is negligible on the scale of the 
fluctuations of E2, and the error in the improved ^(Tr, vr) is essentially the width of the elongated shape in the vertical 
direction. In this case the covariance method leads to error bars « 1/100 of those calculated in the standard way. 

Unfortunate ly, n ot all quantities exhibit a strong covariance with C(1,0). Fig. ^ shows results for the staggered 
susceptibility (p2q) of a 6 x 6 system. In this case there is only a very weak covariance, and hardly any gain in 
accuracy can be achieved. 

It is easy to understand why the covariance with C(l, 0) is particularly strong for ^(Tr, tt) (or indeed any equal-time 
spin correlation): The system is rotationally invariant, but the simulation generates configurations in a representation 
where the z-direction is singled out, and only this component of the correlation function is measured (the other 
components are not easily measurable, which is the case also with standard worldline methods). Measurements based 
on a particular set of configurations (a single bin or a bootstrap sample) will inevitably be affected by some deviations 
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from perfect rotational invariance. This is manifested as amplitude fluctuations in the particular spin component 
measured, and cause the covariance effects seen in the data discussed above. The ability of the local Monte Carlo 
updates to rotate the direction of the antiferromagnetic order in spin space diminishes with increasing size, leading to 
large statistical fluctuations in the conventional estimate of the correlations. The fact that the energy fluctuations do 
not increase with the system size can, in the same way, be traced to the rotationally invariant nature of the estimator 
(\lw. A somewhat more formal discussion of the covariance error reduction scheme can be found in Ref. li'n. 



III. RESULTS 
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Simulations of L x L systems with L < 16 (only even L) were carried out at inverse temperatures /3 = AL and 8L. 
Within statistical errors the results are indistinguishable, indicating that in both cases the ground state completely 
dominates the behavior of the calculated quantities. This can also be checked using the finite-size singlet-triplet gap 
scaling predicted from chiral perturbation theory.E3 For L = 16 and (3 = 128, this gives an estimate of ~ 10""^ for the 
relative error in the calculated ground state energy due to excited states (note that since the simulations are carried 
out in the canonical ensemble, only = states are mixed in). Ear the smaller systems the errors are even smaller 
(the gap scales as 1/A^). All results discussed here are for /3 — 8LE3 

The statistical errors of the calculated energies are as small as ~ 10^^ for all L studied. This accuracy excee^ 
by a factor 5-6 the the most accurate results previously reported for L = 4 — 16; the GFMC results by Runge. 
Comparing the two sets of results, they agree for L < 8, but for the larger sizes the GFMC data is consistently 
higher by 2 — 3 GFMC error bars. Given the agreement to a relative accuracy of less than 10~^ between the SSE 
result for L — 6 and the exact result, and the non- approximate nature of the algorithm, it is hard to see why there 
should be any systematic errors in the SSE data for the larger lattices. Note that any remaining finite-temperature 
effects would lead to an overestimation of the energy, and hence could not explain the discrepancy with the GFMC 
data. As discussed above, care has been taken to verify that in fact the finite temperature effects are well below 
the statistical errors. GFMC calculations, on the other hand, are in general expected to be affected by a small 
systematic error originating from "population control" of the varying number of "random walkers" used in that tatpa 
of simulation. In Runge's calculation, attempts were made to remove such bias more effectively than in previouaiil'EJ 
GFMC calculations. However, a small remaining systematic error could not be ruled out, and the effect was expected 
to be an over-estimation of the energyJlj The discrepancy found here is therefore not completely surprising. The 
energies obtained with the SSE algorithm are listed in Table It is gratifying to note that the SSE energy is indeed 
nicely self-averaging — despite the considerably fewer numbers of MC steps performed for the larger systems the 
relative errors do not differ much from the smaller sizes. m 

Two definitions of the sublattice magnetization have been frequently used in previous studies,™ and will be used 
here as well. The first definition is in terms of the staggered structure factor, 

Afi^(L) = 3S'(7r,7r)/i^ (33) 

and the second one uses the spin-spin correlation function at the largest separation on the finite lattice, 

M|(L) =3C(i/2,i/2). (34) 

The factors 3 are included to account for rotational averaging of the ^-component of the correlations. Mi{L) and 
M2{L) should, of course, scale to the same sublattice magnetization in the thermodynamic limit. Using covariance 
error reduction, both were determined to within statisticals errors of « 10~* (slightly larger for the largest systems). 
This accuracy also exceeds that of previous studies. The results for both ^(Tr, tt) and C(L/2, L/2) are listed i Table Q. 

As discussed in Sec. II, the spin stiffness can be obtained directly bj4 measuring the square of the fluctuating winding 
number. However, as pointed out recently by Einarsson and Schulzp3 the two terms in Eq. ( p6| ) have different leading 
size-corrections; ^ Xjl? for E and ^ 1/L for A^. Therefore, A^ is also calculated separately. There is a small 
discrepancy between the QMC results for L = 4 and L = 6, and the Lanczos results reported in Ref. Adjusting for 
different factors in the definitions, the Lanczos results are As(4) — 0.04832 and As(6) = 0.06723, whereas the QMC 
results obtained here are As(4) = 0.04841(2) and As(6) = 0.06791(3). The reason for the discrepancy is not clear, 
but carrying out 4x4 exact diagonalizations with weak twist-fields included in the Hamiltonian, and subsequently 
calculating the derivative in Eq. (|2C| ) numerically, gives As(4) = 0.04840, in good agreement with the QMC result. 
Hence, there is reason to believe that the QMC results are correct. 

The uniform suscepjtihili|tj4 has typically been obtained in numerical studies via a definition in terms of the singlet- 
triplet excitation gap.E^fEjtj Here a different approach is taken, not requiring simulations in the S = 1 sector. The 
q-dependent susceptibility, 
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N 



is calculated using Eq. (15), and its value at the longest wavelength, qi = 27r/L, is taken as the definition of the 
finite-size uniform susceptibility [due to the finite-size gap and the conserved magnetization, xiQ = 0) of course 
vanishes identically] . In order to give the correct transverse susceptibility of a system with broken symmetry in the 
thermodynamic limit, the result has to be adjusted by factor 3/2. Hence, the definition is 

X±{L) = lx{2n/L,0). (36) 

The spinwave velocity can be obtained from the infinite-size values of ps and x± according to the general hydrodynamic 
relation 



c = VpJxI- (37) 

The above quantities will now be scaled to the thermodynamic limit using g^g to appropriate scaling forms. Chiral 
perturbation theory gives the following scaling behavior for the ground state energy and the sublattice magnetization 
defined according to Eq. (|3^) [parameters without the argument L will henceforth denote the infinite-size values] : 

E{L)^E + (3cl^ + ^^^ + ..., (38a) 
A'P 1 

Mf (L) = M^ + a - + ..., (38b) 

cx± L 



where a = 0.62075 and /? = — 1.4377.r.1 The leading corrections have been obtained also frorcLxenormalization group 
calculations for the nonlinear a model,E§Ej and their orders also agree with spiawaye theory.L^ 

Spinwave theory gives that the leading corrections to and M2 are ~ l/i,E3l3 and this should be the case also 
for x±(^) due to the linear spinwave spectrum for small q. To the author's knowledge, there are no more detailed 
predictions for the scaling behavior of these quantities. 

Individually fitting all the parameters, it is found that the high accuracy of E(L) necessitates the inclusion also of a 



term ^ \/L° in Eq. (38a). Both Mf{L) and Af|(L) require corrections up to order l/L^. The QMC results for As(L) 
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and x±(i) are less accurate, and only linear and quadratic terms are needed. Hence, the following size dependences 
are assumed 

E{L) 
Ml{L) 
MiiL) 

X±{L) 

The predicted scaling forms (^), together with the hydrodynamic relation (^7|) and the expression ( |26| ) for the spin 
stiffness p^, imply the following constraints among the parameters and size-corrections: 

A, = ~{l/3)[E + 2aM^e3/{Pmi)] (40a) 
X± = af3M^/{mie3) (40b) 
64 = mie3/(4a/3Af2). (40c) 

All the scaling forms ( p9[ ) are hence coupled to a high degree, and a good simultaneous fit of all parameters will 
strongly support the field theoretical predictions (|3^). 

Data for all sizes L = 4 — 16 can be included in fits with good values of x^ per degree of freedom (x2/D0F), except 
in the case of X_l(-Zj) for which L ~ 4 has to be excluded (not surprising, since the smallest wave- vector qi used in 
the definition is as large as 7r/2 for L ~ 4). Both As(16) and x_l(16) have error bars too large to be useful, and are 
therefore also excluded. 

Extrapolating the infinite-size parameters from fits to a small number of points, one has to take into account the 
fact that there are higher-order corrections present, which by necessity have been neglected in the scaling forms used. 
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The statistical errors of the extrapolated parameters may be smaller than the systematic errors introduced due to 
this neglect (even though the fit may be good). In order to minimize this type of subtle errors, the L = 4 data were 
excluded from all the fits discussed in the following. This leads to larger statistical fluctuations but should significantly 
reduce the risk of underestimating the errors (the largest neglected correction to is 11 times larger for L — 4 than 
for L — 6, and for the other quantities 3 — 5 times larger). 

Before considering the full fit ( ^9|) with all the constraints (^0|), it is instructive to consider first the results of 
individual, unconstrained fits to all the different quantities. The effects of including the constraints can then be 
judged in light of these results. 

Completely independent fits give E = -0.66943(2), M = 0.3062(6) [from Mi(L)], M = 0.3068(9) [from M2(L)], 
Ps = 0.179(4), and x± = 0.063(1). Note that Mi{L) and M2(L) give the same sublattice magnetization M within 
statistical errors, as they should. Using ( |37| ) the spinwave velocity is c = 1.69(2). These parameters are in goOjd-jgeneral 
agreement with previous calculations, except that the energy is slightly lower than the best GFMC estimatejla due to 
the discrepancies in the firute-sizc data discussed above. The sublattice magnetization is a bit lower than the recent 
result by Beard and WieseN [M = 0.3085(2)1. 

The consistency with the scaling forms (|38| ) can of course also be tested with these independent fits. The leading 
energy correction is found to be 63 = —2.43(11), whereas Eq. (38a) in combination with the above estimate for c 
gives 63 = /3c = 2 .43( 3). The constant of the linear term in Mi{L) is mi — 0.574(9), whereas the right hand side of 
the scali ng fo rm ( |38b| ) gives mi = aM^ / {cx±) = 0.550(10). The subleading energy correction of the fit is 64 — 4(1), 
and Eq. (38a) gives 64 = c^/(4ps) = 4.0(1). Hence, there is good consistency with the predicted scaling forms, within 
a few percent for the leading terms, but with a statistical uncertainty as large as « 25% for the subleading energy 
correction. 

The leading corrections have been derived in several differe nt wa ys.Ejn J'EZl Given then the good numerical agreement 
found above, it is now reasonable to enforce the constraints ([40aD and (40b) which involve these terms. From such a 
fit, with the constraint on the subleading energy correction (40c) left unenforced, one can get a better estimate of the 
size of the subleading term. The constraint that Mi{L) and M2{L) extrapolate to the same M is also enforced. 

One of the early nonlineaK-ff model calculations of the finite-size behavior indicated that the subleading correction to 
E was ^ 1/i^, not ~ 1/L''.c3 Previous numerical calculations were not accurate enough to distinguish between these 
forms. However, Runge noted that the value of c extracted from the leading correction was in rather poor agreement 
with other estimates if a 1/L^ subleading correction was used, and that a slightly better value was obtained using 
Even with the accuracy of the QMC results for E{L) obtained here, individual fits using the two different subleading 
corrections (and including also the next higher-order correction in both cases) cannot by themselves definitely rule 
out the absence of the term, alth ough the fit including it is better. However, it is not at all possible to obtain a 
good fit constrained by (40a) and (40b) without the term, x^/DOF being as high as ~ 50 in this case. With the 
1/L'* term y;^ / DOFk, 0.9. Hence, knowing the constraints on the leading corrections, the present data unambiguously 
require that the subleading energy term is ~ 1/L^. 

The parameters obtained in the partially constrained fit are: E = —0.669436(5), M — 0.3071(3), ps = 0.176(2), 
X± = 0.0623(10), and c = 1.681(14). The statistical errors are here significantly reduced relative to the previous 
unconstrained fits, and the two sets of parameters are consistent with each other. The leading corrections to E and 
M are now, of course, in complete agreement with the theoretical prediction. The subleading energy correction of 
the fit is 64 = 4.17(23), whereas Eq. ( ^8a| ) with the above parameters gives 64 = c^/(4ps) = 4.01(7). Hence, it is 
now also confirmed, at the 5% accuracy level, that the size of the_subleading energy correction agrees with the chiral 
perturbation theory prediction by Hasenfratz and NiedermayerO It is remarkable that the derivation of this very 
detailed theoretical result is based purely on symmetry and dimensionality considerations .E3 

With the subleading energy correction now firmly established, the last constraint (40c) can also be enforced, and 
the parameters of this fit are taken as the final results. This coupled fit still has x^/DOFw 0.9 (the total number 
of parameters is 14, and a total of 28 data points were used). Looking at for the five individual curves, they all 
represent good fits to their respective data sets. Hence, the constrained fit is in all respects a good one. All the ground 
state parameters obtained from the fit are listed in Table III , along with the leading and subleading corrections to the 
energy and the sublattice magnetization. There are only minor changes relative to the previous partially constrained 
fit, the main improvement in-accuracy being for the spinwave velocity. The errors of the parameters were calculated 
using the bootstrap method)^ i.e., fits were carried out for a large number of bootstrap samples of the QMC data, 
and the error is defined as one standard deviation of the parameters of those fits. As explained in Sec. II-B, this 
should be a very accurate method for calculating statistical errors even in highly nonlinear situations such as the 
constrained fit. The QMC data used, along with the fitted curves, are shown in Figs. ^|-^. 

In the figures, it can be observed that even though the L — 4 data are not included in the fit, the fitted curves 
extrapolated to L — 4 are quite close to these QMC points, except in the case of XJ-(4) (for which this point cannot 
even be included in an individual fit). In fact, calculating also the statistical errors of the extrapolations to L = 4 gives 
strong further support to the reliability of the procedures used. For all quantities except the energy, the statistical 
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errors are found to be comparable at L = 4 and L = oo. For the energy the fluctuations are more than 20 times larger 
at L = 4, due to the high order of the leading correction. Both E{A) and As (4) are within one standard deviation of the 
extrapolated results. The sublattice magnetizations Mi (4) and M2(4) both deviate by 2.5 standard deviations, and 
X_L deviates by 3 standard deviations. These rather small deviations clearly indicate that there are only minor effects 
of neglected higher-order corrections. The extrapolations to L = c» should of course be significantly less sensitive 
to systematic errors, since the neglected corrections rapidly vanish for the larger sizes. As already discussed above, 
the neglected corrections are several times larger a.t L — 4 than at L = 6 (the smallest size considered in the fits). 
Hence, any remaining systematic errors in the infinite-size extrapolations listed in Table [II should be well below the 
indicated statistical errors. 



IV. SUMMARY AND DISCUSSION 



An extensive study of the ground state pacameters of the 2D Heisenberg model has been presented. Using the 
Stochastic Series Expansion QMC methocE3'L3 in combination with a data analysis scheme utilizing covariance 
effects,EZI results of unprecedented accuracy were obtained for the ground state energy and the sublattice magne- 
tization for systems of linear dimensions up to L = 16. The long-wavelength susceptibility and the spin stiffness were 
also directly calculated in the simulations. 

The QMG-data was extrapolated to the thermodynamic limit using scaling forms predicted from chiral perturba- 
tion theory,EZI supplemented by higher-order terms necessary to obtain good fits. Both the leading and subleading 
corrections were found to agree in magnitude with the theoretical predictions to within a few percent. This is the 
first numerical verification of the predictions of chiral perturbation theoryEJ to subleading order. 

The ground state energy extracted from the fit is the most accurate estimate obtained to date, and is slightly higher 
than the best Green's functiop-Monte Carlo result £3 This discrepancy is most likely due to a "population control" 
bias in the GFMC calculation.lla The spin stiffness and-tJie sublattice magnetization are both lower than the results of 
a recent low-temperature loop algorithm QMC study.E^ The discrepancy appears to be marginally larger than what 
could be explained by statistical fluctuations alone. For the spinwave velocity the results are consistent with each 
other. J— I 

In the QMC study by Beard and WicscEa the size and temperature dependence of the uniform and staggered 
susceptibilities were fit to scaling forms from chiral perturbation theory. Hence, the underlying theory for analyzing 
the data is the same as used in the present study, but the physical quantities used are different, as is the temperature 
regime (low but finite T versus T = in this study). Since both QMC algorithms are "exact", the discrepancies must 
originate from the scaling procedures. In this paper the effects of neglected higher-order corrections were discussed, 
and attempts were made to minimize these as much as possible. Furthermore, as a quantitative check of remaining 
effects of this nature, the calculated scaling functions were extrapolated to lattices smaller than the smallest size used 
in the fit (L = 6). The close ageement with the actual calculated results, along with the high orders of the largest 
neglected corrections, show that any effects of the higher-order terms on the extrapolations to infinite size should be 
well below the carefully computed statistical errors. It can be noted that Beard and Wiese also included subleading 
corrections Jl3 but not to the same high orders as was necessary in the present study (due to the high accuracy of 
the SSE results for the energy and the sublattice magnetization). Another reason to believe that the present study 
is more reliable is that the fit involves in a direct manner the T — finite-size definitions of the same infinite-size 
parameters sought, not functions of those parameters. _ 

In combination with the covariance error reduction scheme,^!! the SSE method is a very efficient method for cal- 
culating correlation functions of isotropic spin models, as exemplified by the results presented here. The covariance 
scheme is most efficient in cases where there are strong long-ranged correlations (where the "bare" estimator for 
the correlation function does not behave well). It is currently being applied in a study the temperature dependence 
of the corteJation length of the weakly coupled Heisenberg bilayer, for larger lattices and lower temperatures than 
previousljytJ possible. This is motivated by recent results obtained from a mapping to a nonlinear a model for this 
system,Eil predicting a much faster divergence of the correlation length as T — > than for the single layer. The 
SSE algorithm has also proven useful in the case of critical, or near-critical systems, such as the bilayer Heisenberg 
model close to its quantum critical point. More accurate_finite-size scalings than previously reported for this,E3 as 
well as other models exhibiting quantum critical behavior ,Cj are possible with the data enhanced using the covariance 
method. l— . 

QMC methods baaad-pa the loop-cluster algorithm J2j have proven to be very efficient in several studies of = 1/2 
Heisenberg modelsJl3'E3£3 and are clearly more efficient than the SSE method in many cases. For example, it was 
shown here that the covariance method cannot significantly improve the accuracy of calculated static susceptibilities, 
which appear to be very accurately given in loop algorithm simulations .113 The method used here for calculating 
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the spin stiffness directly from the winding number fluctuations would probably also be more accurate with loop 
algorithms. However, it is not clear whether loop algorithms can easily produce more accurate results for equal-time 
correlation functions or energies than those presented in this paper. 

The methods discussed here can also be easily extended to higher-spin models. In fact, the covariance scheme has 
an additional advantage for S > 1/2, in that the on-site correlation (5f)^s known exactly, but fluctuates in the 
simulation and exhibits strong covariance with other correlation functions.Ell Detailed studies of various higher-spin 
models should therefore now be feasible also in 2D. 
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L = 4 L^6 

Bexact -0.701780 -0.678872 

El -0.701777(7) -0.678873(4) 

E'z -0.70177(6) -0.67872(9) 

Soxact 1.47481 2.5180 

51 1.47480(4) 2.51799(6) 

52 1.4747(2) 2.5168(8) 



TABLE I. Comparisons of QMC and exact results for 4x4 and 6x6 lattices. iSexact is the exact result for the ground state 



energy per spin, Ei is the QMC estimate obtained from the average length of the series expansion according to F,< 



Q, and 

E2 is the estimate obtained from the nearest-neighbor correlation function C(1,0) calculated according to Eq. (M). Scxact is 
the exact staggered structure factor, Si is the QMC result with the accur acy i ncreased by using the covariance with E2, and 
S2 is the estimate using directly the sum of the correlation functions, Eq. (32a). The exact L — 6 results are from Ref. El 



L -E S{7v,Tv) C{L/2,L/2) 

4 0.701777(7) 1.47480(4) 0.059872(5) 

6 0.678873(4) 2.51799(6) 0.050856(3) 

8 0.673487(4) 3.7939(2) 0.045867(5) 

10 0.671549(4) 5.3124(3) 0.042851(6) 

12 0.670685(5) 7.0780(7) 0.040873(9) 

14 0.670222(7) 9.090(1) 0.03945(1) 

16 0.669976(7) 11.352(2) 0.03839(2) 



TABLE IL QMC results for the ground state energy, the staggered structure factor, and the spin-spin correlation function at 
distance r = (L/2,L/2), The accuracies of the results for S{n,n) and C{L/2, L/2) were increased by using covariance effects, 
as discussed in Sec. IL The simulations were carried out at inverse temperatures (3 — 8L, and the total number of MC steps 
performed for the different systems were 2.5 x 10** (L = 4), 3 x 10* (i = 6), 10** (-L = 8), 10* (L = 10), 3 x 10^ (L = 12), 
1.5 X 10^ (L = 14), and lO'' (L = 16). 



parameters 


E 


-0.669437(5) 


M 


0.3070(3) 


Ps 


0.175(2) 


X± 


0.0625(9) 


c 


1.673(7) 


size corrections 


63 


-2.405(10) 


64 


4.00(6) 


mi 


0.560(6) 


m2 


1.08(5) 



TABLE in. The ground state parameters and the leading and subleading corrections to E and M, as obtained from a fit 
fully constrained by the predictions of chiral perturbation theory. 
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FIG. 1. The operator sequence truncation vs. the number of Monte Carlo steps performed for a 4 x 4 lattice at inverse 
temperature (3 = 32. The final I after 10^ MC steps was I = 804 (the last increase occurred after 6674 steps). The increment 
used was A; — I/IO. 



FIG. 2. The distribution of the power n of the sampled terms in a 5 x 10 MC step simulation of 4 x 4 lattice at inverse 
temperature /3 — 32, after adjusting I as shown in Fig. 1. The lower histogram is the full distribution. The higher, only partially 
visible histogram is the distribution multiplied by a factor 1000. The cut-off was I — 804, which is significantly larger than the 
largest n sampled. Hence, the truncation has not degraded the accuracy of the simulation. 

FIG. 3. Correlations between the staggered structure factor and the energy as obtained from the nearest-neighbor correlation 
function for a 6 x 6 lattice at /3 = 48. Each point represents a bootstrap sample of QMC bin averages. The solid vertical line 
indicates the exact energy, and the solid horizontal line is the exact structure factor [the result for ^(Tr, tt) is given "only" with 
5 significant digits in Ref. which actually implies a small uncertainty on the scale of this figure]. The vertical dotted and 
dashed lines indicate, respectively, the estimate ± one standard deviation of the energy calculated from the nearest-neighbor 
correlation function and the average of n. The dotted horizontal lines indicate the conventional estimate of the structure factor, 
and the dashed ones the improved result obtained using the covariance method. 

FIG. 4. Correlations between the energy as obtained from the nearest-neighbor correlation function and the staggered 
structure factor for a 16 x 16 lattice at /3 = 128. The vertical line indicates the estimate of the energy from the average of n (the 
error is too small to be seen on this scale). The dotted horizontal lines indicate the value ± one standard deviation of S(7r,7r) 
calculated using the conventional estimator. The dashed, almost indistinguishable lines indicate the improved estimate. 

FIG. 5. Correlations between the staggered susceptibility and the energy as obtained from the nearest-neighbor correlation 
function for a 6 x 6 lattice at j3 — 48. Due to the very weak covariance, only a modest increase in accuracy can be achieved for 
X(7r,7r). 

FIG. 6. The ground state energy vs. on two different scales. The solid circles are the QMC data. On the scale of the 

top graph, the error bars are smaller than half thiwadius of the circles. The smaller circles with error bars in the top graph 
are the GFMC and extrapolated results by RungeJl£l The result of the constrained fit including the L > 4 data is indicated by 
the curves. 

FIG. 7. The sublattice magnetization squared, as defined by Eqs. ( |33| ) and (^^, vs. inverse system size. The solid and open 
circles are the QMC data for Mi (L) and M2 (L) , respectively, and the curves are the results of the constrained fit including the 
L > 4 data. The QMC error bars are much smaller than the circles. 



FIG. 8. The spin current-current correlator, Eq. (^, vs. the inverse system size, along with the result of the coupled fit to 
the L > 4 data. 



FIG. 9. The long-wavelength spin susceptibility [2/3 of x±{L)] vs. inverse system size, along with the result of the constrained 
fit to the L > 4 data. 
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